package com.asa.numanaly;

import java.util.ArrayList;
import java.util.List;

import com.asa.HanShu;
import com.asa.HanShuN;

/**
 * 常微分方程初值问题
 * @author Administrator
 *
 */
public class E {
	
	
	/**
	 * Euler方法 *
	 * 这个正经用得少，但是是很多算法的基础，算是初代原型机
	 * @return
	 */
	public static double[][] chuzhi1(double a,double b,int n,double chuzhitiaojian,HanShu hanshu) {
		double[][] result = new double[n+1][2]; 
		double h = (b-a)/n;
		double t = a;
		double w = chuzhitiaojian;
		
//		System.out.println(t+"  "+w);
		result[0][0] = t;
		result[0][1] = w;
		for (int i = 1; i <= n; i++) {
			
			w = w+ h*hanshu.hanshu(t, w);
			t = a+i*h;
			
//			System.out.println(t+"  "+w);
			result[i][0] = t;
			result[i][1] = w;


		}
		

		return result;

	}
	
	
	/**
	 * Runge-Kutta	4阶龙格库塔 
	 * 数据后几位有些对不上
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @return
	 */
	public static double[][] chuzhi2(double a,double b,int n,double chuzhitiaojian,HanShu hanshu) {
		double[][] result = new double[n+1][2]; 
		double h = (b-a)/n;
		double t = a;
		double w = chuzhitiaojian;
		
//		System.out.println(t+"  "+w);
		result[0][0] = t;
		result[0][1] = w;
		for (int i = 1; i <= n; i++) {


			double k1 = h*hanshu.hanshu(t, w);
			double k2 = h*hanshu.hanshu(t+h/2, w+k1/2);
			double k3 = h*hanshu.hanshu(t+h/2, w+k2/2);
			double k4 = h*hanshu.hanshu(t+h, w+k3);

			w = w+(k1+2*k2+2*k3+k4)/6;
			t = a+i*h;
			result[i][0] = t;
			result[i][1] = w;
			
		}
		

		return result;

	}
	
	
	
	/**
	 * Runge-Kutta-Fehlberg	龙格库塔 ，动态调整步长 *
	 * 
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @return
	 */
	public static double[][] chuzhi3(double a,double b,int n,double chuzhitiaojian,HanShu hanshu,double TOL,double hmax,double hmin) {
		double[][] result = new double[n+1][3]; 
		
		double t = a;
		double w = chuzhitiaojian;
		double h = hmax;
		double FLAG = 1;
		int i = 0;
//		System.out.println(t+"  "+w);
		result[0][0] = t;
		result[0][1] = w;
		while(FLAG==1) {

			double k1 = h*hanshu.hanshu(t, w);
			double k2 = h*hanshu.hanshu(t+h/4, w+k1/4);
			double k3 = h*hanshu.hanshu(t+h*3/8, w+k1*3/32+k2*9/32);
			double k4 = h*hanshu.hanshu(t+h*12/13, w+k1*1932/2197-7200*k2/2197+k3*7296/2197);
			double k5 = h*hanshu.hanshu(t+h, w+k1*439/216-8*k2+k3*3680/513-k4*845/4104);
			double k6 = h*hanshu.hanshu(t+h/2, w-k1*8/27+2*k2-k3*3544/2565+k4*1859/4104-k5*11/40);

			double R = Math.abs(k1/360 - k3*128/4275 -k4*2197/75240 + k5/30 + k6*2/55)/h;
			
			
			if (R<= TOL) {
				
				t = t + h;
				w = w + k1*25/216 + k3*1408/2565 + k4*2197/4104 - k5/5;
				
				result[i][0] = t;
				result[i][1] = w;
				result[i][2] = h;
				
			}
			double o = 0.84*Math.pow(TOL/R, 0.25);
			if (o<=0.1) {
				h = 0.1*h;
			}else if (o>=4) {
				h = 4*h;
			}else {
				h = o*h;
			}
			if (h>=hmax) {
				h = hmax;
			}
			if (t>=b) {
				FLAG = 0;
			}else if (t+h>b) {
				h = b -t;
			}else if (h<hmin) {
				FLAG = 0;
				return null;
			}

		
			i++;
		}
		

		return result;

	}
	
	
	/**
	 * Adams 4阶预测校正算法	*
	 * p264
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @return
	 */
	public static double[][] chuzhi4(double a,double b,int n,double chuzhitiaojian,HanShu hanshu) {
		double[][] result = new double[n+1][2]; 
		double[] ti = new double[4];
		double[] wi = new double[4];

		double h = (b-a)/n;
		ti[0] = a;
		wi[0] = chuzhitiaojian;
		
		result[0][0] = ti[0];
		result[0][1] = wi[0];
		for (int i = 1; i <= 3; i++) {
			double k1 = h*hanshu.hanshu(ti[i-1], wi[i-1]);
			double k2 = h*hanshu.hanshu(ti[i-1]+h/2, wi[i-1]+k1/2);
			double k3 = h*hanshu.hanshu(ti[i-1]+h/2, wi[i-1]+k2/2);
			double k4 = h*hanshu.hanshu(ti[i-1]+h, wi[i-1]+k3);
			wi[i] = wi[i-1]+(k1+2*k2+2*k3+k4)/6;
			ti[i] = a+i*h;
			result[i][0] = ti[i];
			result[i][1] = wi[i];

		}
		for (int i = 4; i <= n; i++) {
			
			double t = a +i*h;
			double w = wi[3]+ h*(55*hanshu.hanshu(ti[3], wi[3])
					-59*hanshu.hanshu(ti[2], wi[2]+37*hanshu.hanshu(ti[1], wi[1]) - 9*hanshu.hanshu(ti[0], wi[0])))/24;
			
			w = wi[3]+ h*(9*hanshu.hanshu(t, w)+19*hanshu.hanshu(ti[3], wi[3]-5*hanshu.hanshu(ti[2], wi[2]) 
					+ hanshu.hanshu(ti[1], wi[1])))/24;

			for (int j = 0; j <= 2; j++) {
				ti[j] = ti[j+1];
				wi[j] = wi[j+1];
			}
			ti[3] = t;
			wi[3] = w;
			
			result[i][0] = t;
			result[i][1] = w;
		}
		

		

		return result;

	}
	
	
	public static double[][] RK4(double h,double[] x,double[] v,HanShu hanshu, int index) {
		double[][] result = new double[2][4]; 
		result[0][0] = v[0];
		result[1][0] = x[0];
		for (int j = 1+index; j <= 3; j++) {


			double k1 = h*hanshu.hanshu(x[j-1], v[j-1]);
			double k2 = h*hanshu.hanshu(x[j-1]+h/2, v[j-1]+k1/2);
			double k3 = h*hanshu.hanshu(x[j-1]+h/2, v[j-1]+k2/2);
			double k4 = h*hanshu.hanshu(x[j-1]+h, v[j-1]+k3);

			v[j] = v[j-1]+(k1+2*k2+2*k3+k4)/6;
			x[j] = x[index]+j*h;
			result[0][j] = v[j];
			result[1][j] = x[j];
			
		}
		

		return result;
	}
	
	
	
	/**
	 * Adams 变步长预测校正算法	*
	 * p269
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @return
	 */
	public static List<double[]> chuzhi5(double a,double b,int n,double chuzhitiaojian,HanShu hanshu,double TOL,double hmax,double hmin) {

		List<double[]> result = new ArrayList<double[]>();
		
		double[] ti = new double[4];
		double[] wi = new double[4];

//		double[] x = new double[4];
//		double[] v = new double[4];

		
		double h = hmax;
		ti[0] = a;
		wi[0] = chuzhitiaojian;

		double FLAG = 1;
		double LAST = 0;

		
		
		
		result.add(new double[]{0.0,ti[0],wi[0],h});
		
		
		RK4(h, ti, wi, hanshu,0);//这个是地址传递
		double NFLAG = 1;//表明是RK4计算的
		int i= 4;
		double t = ti[3]+h;
		
		
		while(FLAG==1) {
			
			double WP = wi[i-1]	+h*(55*hanshu.hanshu(ti[i-1], wi[i-1])
					-59*hanshu.hanshu(ti[i - 2], wi[i- 2]+37*hanshu.hanshu(ti[i-3], wi[i-3]) - 9*hanshu.hanshu(ti[i-4], wi[i-4])))/24;
			
			double WC = wi[i-1]+ h*(9*hanshu.hanshu(t, WP)+19*hanshu.hanshu(ti[i-1], wi[i-1]-5*hanshu.hanshu(ti[i-2], wi[i-2]) 
					+ hanshu.hanshu(ti[i-3], wi[i-3])))/24;
			
			
			
			
			double o = 19*Math.abs(WC - WP)/(270*h);
			
			if (o<=TOL) {//接受结果
				
				wi[i] = WC;
				ti[i] = t;
				
				if (NFLAG == 1) {
					
					for (int j = i-3; j <= i; j++) {
						result.add(new double[]{j,ti[j],wi[j],h});
					}
					
				}else {
					result.add(new double[]{i,ti[i],wi[i],h});

					
				}
				
				
				if (LAST ==1) {
					FLAG = 0;
					t = ti[i-1]+h;
					continue;
				}else {
					
					i = i+1;
					NFLAG = 0;
					if (o<0.1*TOL||ti[i-1]+h>b) {
						double q = Math.pow(TOL/(2*o), .25);
						if (q>4) {
							h = 4* h;
						}else {
							h = q*h;
						}
						if (h>hmax) {
							h = hmax;
						}
						if (ti[i-1]+4*h>b) {
							h = (b-ti[i-1])/4;
							LAST = 1;
						}
						RK4(h, ti, wi, hanshu,i-1);
						NFLAG = 1;
						i = i+3;
					}
					
					
				}
				
				
				
				
				
				
				
			}else {//拒绝结果
				
				double q = Math.pow(TOL/(2*o), .25);
				if (q<0.1) {
					h = 0.1*h;
				}else {
					h = q*h;
				}
				if (h<hmin) {
					FLAG = 0;
				}else {
					if (NFLAG==1) {
						i=i-3;
						RK4(h, ti, wi, hanshu, i-1);
						i = i +3;
						NFLAG = 1;
					}
				}
				t = ti[i-1]+h;
				
			}
			
			
			
			
			
		}
		

		

		return result;

	}
	
	
	
	
	
	
	/**
	 * 外推法 *
	 * p275
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @param TOL
	 * @param hmax
	 * @param hmin
	 * @return
	 */
	public static List<double[]> chuzhi6(double a,double b,int n,double chuzhitiaojian,HanShu hanshu,double TOL,double hmax,double hmin) {
		
		List<double[]> result = new ArrayList<double[]>();

		double[] NK = {0,2,4,6,8,12,16,24,32};
		
		double TO = a;
		double WO = chuzhitiaojian;
		double h = hmax;
		double FLAG = 1;
		
		double[][] Q = new double[7][7];
		
		for (int i = 0; i < 7; i++) {
			for (int j = 0; j <= i; j++) {
				Q[i][j] = (NK[i+1]/NK[j])*(NK[i+1]/NK[j]);
			}
		}
		double[] y = new double[9];

		while(FLAG==1) {
			int k = 1;
			double NFLAG = 0;
			
			while(k<=8 && NFLAG==0) {
				double HK = h/NK[k];
				double T = TO;
				double W2 = WO;
				double W3 = W2+HK*hanshu.hanshu(T, W2);
				T = TO +HK;
				
				for (int j = 1; j <= NK[k]-1; j++) {
					double W1 = W2;
					W2 = W3;
					W3 = W1 + 2*HK*hanshu.hanshu(T,W2);
					T = TO +(j+1)*HK;
					
				}
				
				y[k] = (W3+W2+HK*hanshu.hanshu(T,W3))/2;
				
				if (k>=2) {
					int j = k;
					double v = y[k-1];//这一步的反转太难了，用的时候在完善吧
					while(j>=2) {
						y[j-1] = y[j]+(y[j]-y[j-1])/(Q[k-1][j-1]-1);
						j = j-1;
						
					}
					if (Math.abs(y[1]-v)<=TOL) {//大概是这个意思吧，没验证
						NFLAG =1;
					}	
					
				}
				
				k = k+1;

			}
			
			k = k-1;
			
			
			if (NFLAG ==0) {
				
				h = h/2;
				if (h<hmin) {
					System.out.println("hmin exceeded");
					FLAG = 0;
					
				}
				
				
			}else {
				WO = y[1];
				TO = TO+h;
				
				result.add(new double[]{TO,WO,h});
				
				
				if (TO>= b) {
					FLAG = 0;
				}else if(TO+h>b){
					h = b -TO;
				}else if (k<=3&&h<.5*hmax) {
					h = 2*h;
				}
				
				
				
				
				
			}
			
			
			
			
			
		}
		
		
		
		
		
		return result;
		
		
		
	}

	
	
	
	
	
	
	/**
	 * Runge-Kutta	4阶龙格库塔 
	 * 数据后几位有些对不上
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param hanshu
	 * @return
	 */
	public static double[][] chuzhi7(double a,double b,int n,double[] chuzhitiaojian,HanShuN[] hanshu) {
		
		
		
		double[][] result = new double[n+1][hanshu.length+1]; 
		double h = (b-a)/n;
		double t = a;
		double[] w = new double[hanshu.length];
		result[0][0] = t;
		for (int i = 0; i < w.length; i++) {
			w[i] = chuzhitiaojian[i];
			result[0][i+1] = w[i];
		}
		

		for (int i = 1; i <= n; i++) {
			
			double[][] k = new double[5][hanshu.length];
			
			for (int j = 0; j < hanshu.length; j++) {
				k[1][j] = h*hanshu[j].hanshu(t, w);

			}

			double[] bufferw = new double[w.length];
			for (int j = 0; j < w.length; j++) {
				bufferw[j] = w[j]*k[1][j]/2;
			}
			
			for (int j = 0; j < hanshu.length; j++) {
				k[2][j] = h*hanshu[j].hanshu(t+h/2, bufferw);

			}
			for (int j = 0; j < w.length; j++) {
				bufferw[j] = w[j]*k[2][j]/2;
			}
			for (int j = 0; j < hanshu.length; j++) {
				k[3][j] = h*hanshu[j].hanshu(t+h/2, bufferw);

			}
			for (int j = 0; j < w.length; j++) {
				bufferw[j] = w[j]*k[3][j];
			}
			
			for (int j = 0; j < hanshu.length; j++) {
				k[4][j] = h*hanshu[j].hanshu(t+h, bufferw);

			}


			for (int j = 0; j < bufferw.length; j++) {
				w[j] = w[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6;

			}
			
			t = a+i*h;
			double[] r = new double[hanshu.length+1];
			r[0] = t;
			for (int j = 0; j < hanshu.length; j++) {
				r[j+1] = w[j];
			}
			result[i] = r;
		}
		

		return result;

	}
	
	
	
	
	/**
	 * 具有Newton 迭代的梯形法
	 * @param a
	 * @param b
	 * @param n
	 * @param chuzhitiaojian
	 * @param TOL
	 * @param hanshu
	 * @param m
	 * @return
	 */
	public static double[][] chuzhi8(double a,double b,int n,double chuzhitiaojian,double TOL,HanShu hanshu,int m) {
		double[][] result = new double[n+1][2]; 
		double h = (b-a)/n;
		double t = a;
		double w = chuzhitiaojian;
		
		
//		System.out.println(t+"  "+w);
		result[0][0] = t;
		result[0][1] = w;
		for (int i = 1; i <= n; i++) {


			double k1 = w + h*hanshu.hanshu(t, w)/2;
		
			double w0 = k1;
			int j = 1;
			double FLAG = 0;
			
			while(FLAG==0) {
				
				w = w0 -(w0-h*hanshu.hanshu(t+h, w)/2- k1)/(1 - h*hanshu.hanshu(t+h,w0)/2);
				
				if (Math.abs(w-w0)<TOL) {
					FLAG = 1;
				}else {
					j = j+1;
					w0 = w;
					if (j>m) {
						System.out.println("the maximum number of iterations exceeded");
						return null;
					}
				}
				
			}
			
			t = a+i*h;
			result[i][0] = t;
			result[i][1] = w;
		}
		

		return result;

	}
	
	
	
	
	
	
	
	
	
	
	
	
	public static void main(String[] args) {
		
		
		
		HanShu hanshu = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				
				return x - t*t+1;
			}
			
			@Override
			public double hanshu(double x0) {
				// TODO Auto-generated method stub
				return 0;
			}
		};
		double[][] chuzhi1 = chuzhi1(0, 2, 10, 0.5, hanshu );
		
		for (int i = 0; i < chuzhi1.length; i++) {
			for (int j = 0; j < chuzhi1[i].length; j++) {
				System.out.print(chuzhi1[i][j]+"   ");
			}
			System.out.println();
		}
		System.out.println("==============================================");
		double[][] chuzhi2 = chuzhi2(0, 2, 10, 0.5, hanshu );

		for (int i = 0; i < chuzhi2.length; i++) {
			for (int j = 0; j < chuzhi2[i].length; j++) {
				System.out.print(chuzhi2[i][j]+"   ");
			}
			System.out.println();
		}
		
		
		
		double[] x = new double[4];
		
		for (int i = 0; i < x.length; i++) {
			System.out.println(x[i]);
		}
		test(x);
		for (int i = 0; i < x.length; i++) {
			System.out.println(x[i]);
		}
		
		
	}
	
	public static void test(double[] x) {
		
		for (int i = 0; i < x.length; i++) {
			x[i] = 1;
		}
	}
	
	
	
	
	

}
